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Abstract. We describe our method to construct line blanketed NLTE 
model atmospheres for hot stars. We employ the Accelerated Lambda 
Iteration and use statistical methods to deal with metal line blanketing. 



1. Introduction 

Stellar atmospheres are open systems and thus cannot be in thermodynamic 
equilibrium (TE). The "Local Thermodynamic Equilibrium" (LTE) is a work- 
ing hypothesis which assumes TE not for the atmosphere as a whole but for small 
volume elements. As a consequence, the atomic population numbers are depend- 
ing only on the local (electron) temperature and electron density via the Saha- 
Boltzmann equations. Computing models by replacing the Saha-Boltzmann 
equations by the rate equations (statistical equilibrium) are called non-LTE (or 
NLTE) models. NLTE calculations are more costly than LTE calculations, how- 
ever, it is hard to predict if NLTE effects are important in a specific problem. 
Generally, NLTE effects are large at high temperatures and low densities, which 
implies intense radiation fields hence frequent radiative processes and less fre- 
quent particle collisions which tend to enforce LTE conditions. 

We will restrict ourselves here to the classical model atmosphere problem, 

1. e. the solution of the radiation transfer equations assuming hydrostatic, ra- 
diative and statistical equilibrium. The numerical problem going from LTE to 
realistic NLTE models has been solved only recently and is the topic of this 
paper. This was achieved by the development of new numerical techniques for 
model construction and on the availability of atomic data for many species. The 
replacement of the Saha-Boltzmann equations by the atomic rate equations re- 
quires a different numerical solution technique, otherwise metal opacities cannot 
be accounted for at all. Such techniques were developed with big success during 
the last decade, triggered by important papers by Cannon (1973) and Scharmer 
(1981). The Accelerated Lambda Iteration (ALI) is at the heart of of this de- 
velopment. Combined with statistical methods we are finally able to compute 
metal line blanketed NLTE models with a very high level of sophistication. 

2. Overview of problem and solution method 

We assume plane-parallel geometry, which is well justified for most stars because 
the atmospheres are thin compared to the stellar radius. The only parameters 
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which characterize uniquely such an atmosphere are the effective temperature 
(T e ff), which is a measure for the amount of energy transported through the 
atmosphere per unit area and time, the surface gravity (g), and the chemical 
composition. To construct model atmospheres we have to solve simultaneously a 
set of equations that is highly coupled and non-linear. Because of the coupling, 
no equation is determining uniquely a single quantity - all equations determine 
a number of state parameters. However, each of them is usually thought of as 
determining a particular quantity. These equations are: 

• The radiation transfer equations which are solved for the (angular) mean 
intensities Ji,i = 1, . . . , NF, on a pre-chosen frequency grid comprising 
NF points. The formal solution is given by J = AS, where S is the source 
function (Eq.[To|). Although A is written as an operator, one may think of 
A as a process of obtaining the mean intensity from the source function. 

• The hydrostatic equation which determines the total particle density N. 

• The radiative equilibrium equation from which the temperature T follows. 

• The particle conservation equation, determining the electron density n e . 

• The statistical equilibrium equations which are solved for the population 
densities rii,i = 1, . . . , NL of the atomic NLTE levels. 

• The definition equation for a fictitious massive particle density njj which 
is introduced for a convenient representation of the solution procedure. 

This set of equations has to be solved at each point d of a grid comprising ND 
depth points. Thus we are looking for solution vectors 

ip' d = (ni,...,n NL ,n e ,T,n H ,N,J 1 ,...,J NF ), d=l,...,ND. (1) 

The Complete Linearization (CL) method (Auer & Mihalas 1969) solves this 
set by linearizing the equations with respect to all variables. The basic advan- 
tage of the ALI (or "operator splitting") method is that it allows to eliminate 
at the outset the explicit occurrence of the mean intensities J, from the solu- 
tion scheme by expressing these variables by the current, yet to be determined, 
occupation densities and temperature. This is accomplished by an iteration 
procedure which may be written as (suppressing indices indicating depth and 
frequency dependency of variables): 

J n = A*5 n + (A-A*)5 n_1 . (2) 

This means that the actual mean intensity at any iteration step n is computed 
by applying an approximate lambda operator (ALO) A* on the actual (thermal) 
source function S n plus a correction term that is computed from quantities 
known from the previous iteration step. This correction term includes the exact 
lambda operator A which guarantees the exact solution of the radiation transfer 
problem in the limit of convergence: J = AS. The use of A in Eq.|2| only 
indicates that a formal solution of the transfer equation is performed but in our 
application the operator is not constructed explicitly. Instead we employ the 
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Feautrier solution scheme (Mihalas 1978) or a short characteristic method (Olson 
& Kunasz 1987) to solve the transfer equation that is set up as a differential 
equation. The resulting set of equations for the reduced solution vectors 

ij> d = (ni 1 ...,n NL ,n e ,T,nH,N), d=l,...,ND (3) 

is of course still non-linear. The solution is obtained by linearization and itera- 
tion which is performed either with a usual Newton-Raphson iteration or other 
methods (Sect. 5.2.). First model atmosphere calculations with the ALI method 



were performed by Werner (1986). 

Another advantage of the ALI method is that explicit depth coupling of 
the solution vectors Eq.[l] through the transfer equation can be avoided if one 
restricts to diagonal (i.e. local) ALOs. Then the solution vectors Eq.[| are inde- 
pendent from each other and the solution procedure within one iteration step of 
Eq.|2]is much more straightforward. Depth coupling is provided by the correction 
term that involves the exact solution of the transfer equation. The hydrostatic 
equation which also gives an explicit depth coupling, may be taken out of the set 
of equations and can - as experience shows - be solved in between two iteration 
steps of Eq. ^. Then full advantage of a local ALO can be taken. The linearized 
system may be written as: 

iP d = iP° d + 5^ d , (4) 

where ij) d is the current estimate for the solution vector at depth d and 6ij) d is the 
correction vector to be computed. Using a tri-diagonal A* operator the resulting 
system for 5if) d is - like in the classical CL scheme - of block tri-diagonal form 
coupling each depth point d to its nearest neighbors d ± 1 : 

ld^d-i + PdWd + a d5tp d+1 = c d . (5) 

The quantities a, j3, y are (AW x AW) matrices where NN is the total number 
of physical variables, i.e., NN = NL + 4, and c d is the residual error in the 
equations. The solution is obtained by the standard Feautrier scheme. As 
already mentioned, the system Eq. [| breaks into ND independent equations 
5ip d = /3 d 1 c d (d = 1, . . . , ND) when a local A* operator is used. The additional 
numerical effort to set up the subdiagonal matrices and matrix multiplications 
in the tri-diagonal case is outweighed by the faster global convergence of the ALI 
cycle, accomplished by the explicit depth coupling in the linearization procedure 
(Werner 1989). 

The principal advantage of the ALI over the CL method becomes clear at 
this point. Each matrix inversion necessary to solve Eq.[5] requires (NL + 4) 3 
operations whereas in the CL method (NL + NF + 4) 3 operations are needed. 
Since the number of frequency points NF is much larger than the number of 
levels NL, the matrix inversion in the CL approach is dominated by NF. 

Recent developments concern the problem that the total number of atomic 
levels tractable in NLTE with the ALI method described so far is restricted to 
the order of 250, from experience with our model atmosphere code PR02. This 
limit is a consequence of the non-linearity of the equations, and in order to 
overcome it, measures must be taken in order to achieve a linear system whose 
numerical solution is much more stable. Such a pre-conditioning procedure has 
been first applied in the ALI context by Werner & Husfeld (1985) using the core 
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saturation method (Rybicki 1972). More advanced work achieves linearity by 
replacing the A operator with the ^ operator (and by judiciously considering 
some populations as "old" and some as "new" ones within an ALI step) which 
is formally defined by writing: 

Ju = ^vVv , i-e. IV = K/Xu , (6) 



where the total opacity \ v (as defined in Sect, p. 6. ) is calculated from the pre- 



vious ALI cycle. The advantage is that the emissivity r] u (Sect. |3^6] ) is linear in 
the populations, whereas the source function S v is not. Hence the new opera- 
tor ^ gives the solution of the transfer problem by acting on a linear function. 
This idea is based on Rybicki & Hummer (1991) who applied it to the line 
formation problem, i.e. restricting the set of equations to the transfer and rate 
equations and regarding the atmospheric structure as fixed. Hauschildt (1993) 
and Hauschildt & Baron (1999) generalized it to solve the full model atmosphere 
problem. In addition, splitting the set of statistical equations and solving it sep- 
arately for each chemical element means that now many hundreds of levels per 
species are tractable in NLTE. 



3. Basic equations 

We give a short overview of the basic equations to be solved. More detailed 
presentations may be found in Werner & Dreizler (1999) or Werner et al. (2003). 

3.1. Radiation transfer 

Any numerical method requires a formal solution (i.e. atmospheric structure 
already given) of the radiation transfer problem. The radiation transfer at any 
particular depth point can be described by the following equation, formally 
written for positive and negative [i (which is the cosine of the angle between 
direction of propagation and outward directed normal to the surface) separately, 
i.e. for inward and for outward directional intensities I with frequency v: 

± V dIV frf ) =S V -I„(±ri, M e[0,l]. (7) 

t u is the optical depth (which can be defined via the column mass m that is 
used in the other structural equations by dr u = dmxu/p, with the mass density 
p) and S v is the local source function. Introducing the Feautrier intensity 

u u ^ = {I u {p) + I u {-^)) /2 (8) 

we obtain the second-order form (Mihalas 1978): 

P 2 -^ = ^-S„ ne[0,l]. (9) 

We may separate the Thomson emissivity term (scattering from free electrons, 
assumed coherent, with cross-section a e ) from the source function so that 

S u = S' u + n e a e J v /xv, (10) 
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where S' v is the ratio of thermal emissivity to total opacity (Sect. |3.6. ): S' v = 
f}u/Xv Since the mean intensity is the angular integral over the Feautrier inten- 
sity the transfer equation becomes 



2 d^u u ^ 



s' 



n P o~ P 



Xv 



u. 



d\i. 



(11) 



Thomson scattering complicates the situation by explicit angle coupling. As- 
suming complete frequency redistribution in spectral lines, no explicit frequency 
coupling occurs so that the parallel solution for all frequencies enables a very 
efficient vectorization on the computer. 

3.2. Statistical equilibrium 

Rate equations For each atomic level i the rate equation describes the equilib- 
rium of rates into and rates out of this level: 



n 



E^-E^ = °- 



(12) 



The rate coefficients Pij have radiative and collisional components: Pij = Rij + 
Cij. The radiative downward rate for example is given by: 



R 



Hi 



4tt 



hv 



+ J v \e- hv l kT dv. 



(13) 



Photon cross-sections are denoted by aij(y). (rii/rij)* is the Boltzmann LTE 
population ratio. The computation of collisional rates is generally dependent on 
the specific ion or even transition. 

Abundance definition equation The rate equation for the highest level of a 
given chemical species is redundant. It is replaced by the abundance definition 
equation. Summation over all levels usually includes not only NLTE levels but 
also levels which are treated in LTE, according to the specification in the model 
atom. Denoting the number of ionization stages of species k with NION(k), the 
number of NLTE and LTE levels per ion with NL(l) and LTE (I), respectively, 
and the number fraction of species k with y^, we can write: 



NION(k) 

E 

i=i 



NL(l) LTE(l) 
E n kU+ E 



kit 



=1 



=1 



y k (N - n e 



(14) 



Charge conservation We close the system of statistical equations by invok- 
ing charge conservation. We denote the total number of chemical species with 
NATOM, the charge of ion I with q(l) (in units of the electron charge) and write: 



NATOM NION(k) 

E E 9(0 



k=l 



1=1 



NL(l) LTE (I) 

v " n ku+ E n iu 



8=1 



(15) 
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Complete statistical equilibrium equations We introduce a vector comprising 
the occupation numbers of all NLTE levels, n = (m, . . . , Wjvx)- Then the sta- 
tistical equilibrium equation is written as: 

An = b. (16) 

The gross structure of the rate matrix A is of block matrix form, because tran- 
sitions between levels occur within one ionization stage or to the ground state 
of the next ion. The structure is complicated by ionizations into excited levels 
and by the abundance definition and charge conservation equations which give 
additional non-zero elements in the corresponding lines of A. 

3.3. Radiative equilibrium 

Radiative equilibrium can be enforced by adjusting the temperature stratifica- 
tion either during the linearization procedure or in between ALI iterations. In 
the former case a linear combination of two different formulations can be used 
and in the latter case the classical Unsold-Lucy temperature correction proce- 
dure (Lucy 1964) is utilized. The latter is particularly interesting, because it 
allows to exploit the blocked form of the rate coefficient matrix. This will enable 
an economic block-by-block solution followed by a subsequent Unsold-Lucy tem- 
perature correction step. On the other side, however, this correction procedure 
may decelerate the global convergence behavior of the ALI iteration. 

The two forms of expressing the radiative equilibrium condition follow from 
the postulation that the energy emitted by a volume element per unit time is 
equal to the absorbed energy per unit time (integral form): 



oc 







Xu(S v - J„) du = 0. (17) 



This formulation is equivalent to invoking flux constancy throughout the atmo- 
sphere (differential form) involving the nominal flux TL: 

roc Q 

/ 7 r-(fvJ„)dv-H = 0, (18) 

JO OT v 

where /„ is the variable Eddington factor, defined as 

fu= [ V 2 u V u,dn/ [ u u ^dfj, (19) 
Jo 1 Jo 

and computed from the Feautrier intensity u Vfl (Eq.|8|) after the formal solution. 
The differential form is more accurate at large depths, while the integral form 
behaves numerically better at small depths. Instead of arbitrarily selecting that 
depth in the atmosphere where we switch from one formulation to the other, 
we use a linear combination of both constraint equations which guarantees a 
smooth transition with depth, based on physical grounds: 



— / Xu{S v - Ju) dv + A*j T—ihJu) dv - A*jH -F = 0. (20) 

Kj Jo JO OTv 



For details on this equation and on our implementation of the Unsold-Lucy 
procedure see Dreizler (2003). 



NLTE Radiative Transfer and Model Atmospheres of Hot Stars 



7 



3.4. Hydrostatic equilibrium 

The total pressure P comprises gas, radiation and turbulent pressures, so that 
we write the equation for hydrostatic equilibrium as: 



±P=± 
dm dm 



(Air f 00 1 \ 

NkT f„J v dv + -pvl rh j = g (21) 



with Boltzmann's constant k and the turbulent velocity Uturb- The hydrostatic 
equation may either be solved simultaneously with all other equations or sep- 
arately in between iterations. The overall convergence behavior is usually the 
same in both cases. If taken into the linearization scheme and a local A* oper- 
ator is used then, like in the case of the radiative equilibrium equation, explicit 
depth coupling enters via the depth derivative d/dm. Again, solution of the 
linearized equations has to proceed inwards starting at the outer boundary. 

3.5. Particle conservation and fictitious massive particle density 

The total particle density N is the sum of electron density plus the population 
density of all atomic states, LTE and NLTE levels: 



NATOM NION(k) 

N = n e+ £ £ 
k=l 1=1 



NL(l) LTE (I) 

Y, n kli + E n kU 
i=l i=l 



(22) 



A fictitious massive particle density nn is defined by: 



NATOM NION(k) 

n H = E m k E 



k=l 



1=1 



NL(l) LTE(l) 

E n kli+ E n kli 
i=l i=l 



(23) 



The mass of a chemical species in AMU is denoted by m k . Introducing the mass 
of a hydrogen atom mn, we then may simply write for the material density 



p = n H m H - 



(24) 



3.6. Opacity and emissivity 

Thermal opacity and emissivity are made up by radiative bound-bound, bound- 
free and free-free transitions. For each species we compute and sum up: 



NL{1) NL(l) 



NION 

E 

1=1 |_ i=1 3>* 

NL(l) NL(l+l) 

+ E E ^K-J+L 
i=l j>i 



-hiy-Vi^jkT 



9lj 
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/NL(l+l) 
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LTE(l+l) 
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(25) 
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where the total opacity includes Thomson scattering, i.e. Xv = K u + n e cr e , and 
1)v 



2huy c 2 



NION 

E 



1=1 



NL(l) NL(l) 

E E vu^n^e-^-^ (26) 

i=l j>i 9l 3 



NL(l) NLQ+l) 
+ E E ^i + i, k (v)<e- h ^ kT 

i=l j>i 

(NL(l+l) LTE(l+l) 
E n i+i,i+ E n m, 
i=l i=l 

CU-+1+1 k( u ) denotes the cross-section for photoionization from level i of ion I into 
level k of ion I + 1. The double summation over the bound- free continua takes 
into account the possibility that a particular level may be ionized into more than 
one level of the next high ion. The source function used for the approximate 
radiation transfer is r) v /n v , thus, excludes Thomson scattering. For the exact 



formal solution of course, the total opacity Xv in the expression Eq. 1C includes 
the Thomson term (n e a e ). 

4. The Accelerated Lambda Iteration (ALI) 

In all constraint equations described above the mean intensities J u are sub- 
stituted by the approximate radiation field Eq.[2] in order to eliminate these 
variables from the solution vector Eq. |]]. In principle the ALO may be of ar- 
bitrary form as long as the iteration procedure converges. In practice however 
an optimum choice is desired in order to achieve convergence with a minimum 
amount of iteration steps. 

In the case of diagonal (local) lambda operators the mean intensity J d at a 
particular depth d in the current iteration step is computed solely from the local 
source function and a correction term AJ^, the latter involving the source 
functions (of all depths) from the previous iteration. Dropping the iteration 
count and introducing indices denoting depth points we can rewrite Eq.|2|: 

J d = A\ d S d + AJ d . (27) 

In the discrete form we now think of A* as a matrix acting on a vector whose 
elements comprise the source functions of all depths. Then A^ d is the diagonal 
element of the A* matrix corresponding to depth point d. Writing A dd = B d we 
have a purely local expression for the mean intensity: 

Jd = B d S d + AJ d . (28) 

Much better convergence is obtained if the mean intensity is computed not only 
from the local source function but also from the source function of the neighbor- 
ing depths points. Then the matrix representation of A* is of tri-diagonal form 
and we may write: 



Jd = Cd-iSd-i + B d S d + A d+ iS d+ i + A J d , 



(29) 
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where Cd-i and ^4^+1 represent the upper and lower subdiagonal elements of A* 
and Sd±i the source functions at the adjacent depths. We emphasize again that 
the actual source functions in Eq. ^ are computed from the actual population 
densities and temperature which are unknown. We therefore have a non-linear 
set of equations which is solved by either a Newton-Raphson iteration or other 
techniques, resulting in the solution of a tri-diagonal linear equation of the form 
Eq.§ It was shown in Olson et al. (1986) that the elements of the optimum A* 
matrix are given by the corresponding elements of the exact A matrix. 



5. Solution of the non-linear equations by iteration 

The complete set of non-linear equations for a single iteration step Eq. [| com- 
prises at each depth the equations for statistical, radiative, and hydrostatic 
equilibrium and the particle conservation equation. For the numerical solution 
we introduce discrete depth and frequency grids. The equations are then lin- 
earized and solved by a suitable iterative scheme. Explicit angle dependency of 
the radiation field is not required here and consequently eliminated by the use of 
variable Eddington factors. Angle dependency is only considered in the formal 
solution of the transfer equation. 



5.1. Discretization and Linearization 

A depth grid is set up and we start from a gray approximation computing a 
LTE continuum model using the Unsold-Lucy temperature correction proce- 
dure. Depth points (typical number is 90) are set equidistantly on a logarithmic 
(Rosseland) optical depth scale. Our NLTE code uses the column mass as an 
independent depth variable. 

The frequency grid is established based upon the atomic data input file. 
Frequency points are set blue- and redward of each absorption edge and for 
each spectral line. Gaps are filled in by setting continuum points. Finally, 
the quadrature weights are computed. Frequency integrals appearing e.g. in 



Eq. 2C are replaced by quadrature sums and differential quotients involving depth 
derivatives by difference quotients. 

All variables x are replaced by x — > x + 5x where 8x denotes a small per- 
turbation of x. Terms not linear in these perturbations are neglected. The 
perturbations are expressed by perturbations of the basic variables: 

Q O O O NL r\ 

OX rm OX . OX . llT OX . OT r . . 

ix ~ af ST + a^ Sn < + m iN + a^ Sn » + £ ^ ,n >- (30> 



5.2. Newton-Raphson iteration and alternative methods 

As described in Sect. |2] the linearized equations have a tri-diagonal block-matrix 
form, see Eq.||. Inversion of the grand matrix (= T sized (NN ■ ND) x (NN ■ 
ND), i.e. about 10 4 x 10 4 in typical applications) is performed with a block- 
Gaussian elimination scheme, which means that our iteration of the non-linear 
equations represents a multi-dimensional Newton-Raphson method. The prob- 
lem is structurally simplified when explicit depth coupling is avoided by the use 
of a local ALO, however, the numerical effort is not much reduced, because in 
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both cases the main effort hes with the inversion of matrices sized NN x NN. 
The Newton-Raphson iteration involves two numerically expensive steps, first 
setting up the Jacobian (comprising a, (3, 7) and then inverting it. Additionally, 
the matrix inversions necessary to solve Eq. |g limit their size to about NN = 250 
because otherwise numerical accuracy is lost. 

Two variants recently introduced in stellar atmosphere calculations are able 
to improve both, numerical accuracy and, most of all, computational speed. 
Broyden's (1965) variant belongs to the family of quasi-Newton methods and it 
was first used in model atmosphere calculations in Dreizler & Werner (1991), 
Hamann etal. (1991) and Koesterke etal. (1992). It avoids the repeated set-up 
of the Jacobian by the use of an update formula. On top of this, it also gives 
an update formula for the inverse Jacobian. Another variant, the Kantorovich 
method was introduced into model atmosphere calculations by Hubeny & Lanz 
(1992). This method simply keeps fixed the Jacobian during the linearization 
cycle and it is surprisingly stable. 



6. NLTE metal line blanketing 

Despite the capacity increase for the NLTE treatment of model atmosphere prob- 
lems by introducing the ALI method combined with pre-conditioning techniques, 
the blanketing by millions of lines from the iron group elements arising from tran- 
sitions between some 10 5 levels could only be attacked with the help of statistical 
methods. These have been introduced into NLTE model atmosphere work by 
Anderson (1989). At the outset, model atoms are constructed by combining 
many thousands of levels into a relatively small number of superlevels which can 
be treated with ALI (or other) methods. Then, in order to reduce the computa- 
tional effort, two approaches were developed which vastly decrease the number 
of frequency points (and hence the number of transfer equations to be solved) to 
describe properly the complex frequency dependence of the opacity. These two 
approaches have their roots in LTE modeling techniques, where for the same 
reason statistical methods are applied for the opacity treatment: The Opacity 
Distribution Function (ODF) and Opacity Sampling (OS) approaches. Both are 
based on the circumstance that the opacity (in the LTE approximation) is a 
function of two only local thermodynamic quantities. Roughly speaking, each 
opacity source can be written in terms of a population density and a photon 
cross-section for the respective radiative transition: k u ~ niai u (v). In LTE the 
population follows from the Saha-Boltzmann equations, hence ri[ = ni(n e ,T). 
The OS and ODF methods use such pre-tabulated (on a very fine frequency 
mesh) K u (n e ,T) during the model atmosphere calculations. The NLTE situa- 
tion is more complicated, because pre-tabulation of opacities is not useful. The 
population densities at any depth now also depend explicitly on the radiation 
field (via the rate equations which substitute the TE Saha-Boltzmann statis- 
tics) and thus on the populations in each other depth of the atmosphere. As a 
consequence, the OS and ODF methods are not applied to opacity tabulations, 
but on tabulations of the photon cross-sections These do depend on local 

quantities only, e.g. line broadening by Stark and Doppler effects is calculated 
from T and n e . In the NLTE case the cross-section takes over the role which 
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the opacity played in the LTE case. More details on this topic can be found in 
Werner k Dreizler (1999). 

Our strategy is the following. Before any model atmosphere calculation 
is started, the atomic data are prepared by constructing super levels, and the 
cross-sections for superlines. Then these cross-sections are either sampled on a 
coarse frequency grid or ODFs are constructed. These data are put into the 
model atom which is read by the model code. 

7. Summary 

The construction of metal line blanketed models in hydrostatic and radiative 
equilibrium under NLTE conditions was the last and long-standing problem of 
classical model atmosphere theory and it is finally solved with a high degree of 
sophistication. The essential milestones for this development, starting from the 
pioneering work of Auer & Mihalas (1969) are: 

• Introduction of the Accelerated Lambda Iteration (ALI, or "operator split- 
ting" methods), based upon early work by Cannon (1973) and Scharmer 
(1981). First ALI model atmospheres were constructed by Werner (1986). 

• Introduction of statistical approaches to treat the iron group elements in 
NLTE by Anderson (1989). 

• Linear formulation of the statistical equations (Rybicki & Hummer 1991, 
Hauschildt 1993). 

• Computation of atomic data by Kurucz (1991), by the Opacity Project 
(Seaton etal. 1994) and the Iron Project (Hummer etal. 1993). 
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